Méthode 1 :
soit $X\sim\mathcal{G}(p)$ qui suit une loi géométrique de paramètre 𝑝.
$$ \forall k \in \mathbb{N}^*, \ \ \ P(X = k) = (p-1)^{k-1}\cdot p = q^{k-1}\cdot p $$C'est la loi d'attente d'un premier succès dans une suite d'épreuves répétées indépendantes ayant chacun la probabilité de succès 𝑝
Par conséquent si : $$ (U_i)_{i\geq 1} \sim^{i.i.d} \mathcal{U}(0,1) $$
Alors la variable : $$ Y := \min\{k \in \mathbb{N}^* : U_k < p\} \sim \mathcal{G}(p) $$
Méthode 2 :
code :
proposition de générateur de la loi géométrique
import numpy as np
# méthode 1
def rand_geom(n, p):
"""
n : taille de l'échantillon
p : paramètre de la loi géométrique
"""
rg = np.zeros(n)
for i in range(n):
ru = np.random.rand()
k = 1
while ru > p:
k += 1
ru = np.random.rand()
rg[i] = k
return rg
# méthode 2
def rand_geom2(n, p):
"""
n : taille de l'échantillon
p : paramètre de la loi géométrique
"""
ru = np.random.rand(n)
X = np.log(ru)/np.log(1-p)
Y = np.floor(X)+1
return(Y)
def proba_geom(k, p):
return((1-p)**(k-1) * p)
générons un échantillon de taille (exemple) $n = 10000$ avec le paramètre $p = 0.3$
n = 10000
p = 0.3
sample_geom = rand_geom(n, p)
sample_geom2 = rand_geom2(n, p)
affichons son histogramme ainsi que la vraie distribution de la loi géométrique
import matplotlib.pyplot as plt
plt.hist(sample_geom, color = "steelblue",
density = True,
bins = int(sample_geom.max()),
alpha = 0.5,
label = "Notre échantillon (1)")
plt.hist(sample_geom2, color = "red",
density = True,
bins = int(sample_geom2.max()),
alpha = 0.5,
label = "Notre échantillon (2)")
x = np.arange(1, max(sample_geom.max(), sample_geom2.max()))
y = proba_geom(x, p)
plt.plot(x, y, color = "orange", label = "Vrai proba", marker = "o")
plt.xlabel("k")
plt.ylabel("probabilité")
plt.title(f"histogramme et pdf de la loi géométrique : n = {n} et p = {p}")
plt.legend()
plt.show()
Regardons les deux extrêmes : $p = 0.1$ et $p=0.9$
n = 10000
p = 0.1
sample_geom = rand_geom(n, p)
sample_geom2 = rand_geom2(n, p)
plt.hist(sample_geom, color = "steelblue",
density = True,
bins = int(sample_geom.max()),
alpha = 0.5,
label = "Notre échantillon (1)")
plt.hist(sample_geom2, color = "red",
density = True,
bins = int(sample_geom2.max()),
alpha = 0.5,
label = "Notre échantillon (2)")
x = np.arange(1, max(sample_geom.max(), sample_geom2.max()))
y = proba_geom(x, p)
plt.plot(x, y, color = "orange", label = "Vrai proba", marker = "o")
plt.xlabel("k")
plt.ylabel("probabilité")
plt.title(f"histogramme et pdf de la loi géométrique : n = {n} et p = {p}")
plt.legend()
plt.show()
n = 10000
p = 0.9
sample_geom = rand_geom(n, p)
sample_geom2 = rand_geom2(n, p)
plt.hist(sample_geom, color = "steelblue",
density = True,
bins = int(sample_geom.max()),
alpha = 0.5,
label = "Notre échantillon (1)")
plt.hist(sample_geom2, color = "red",
density = True,
bins = int(sample_geom2.max()),
alpha = 0.5,
label = "Notre échantillon (2)")
x = np.arange(1, max(sample_geom.max(), sample_geom2.max()))
y = proba_geom(x, p)
plt.plot(x, y, color = "orange", label = "Vrai proba", marker = "o")
plt.xlabel("k")
plt.ylabel("probabilité")
plt.title(f"histogramme et pdf de la loi géométrique : n = {n} et p = {p}")
plt.legend()
plt.show()
code
Simulons notre vecteur gaussien de dimension $d\geq2$. Nous pouvons utiliser l'exercice 3 pour les simuler mais il est aussi possible de les simuler via grâce à la librairie numy : numpy.random.multivariate_normal
import numpy as np
d = 3
n = 5000
mu = np.zeros(d)
sigma = np.eye(d)
rn = np.random.multivariate_normal(mu, sigma, n)
Jettons un coup d'oeil à l'échantillon généré (vous pouvez utiliser le module mplot3d de matplotlib)
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = Axes3D(fig)
ax.scatter(rn[:,0], rn[:,1], rn[:,2])
ax.set_xlabel('$x_1$')
ax.set_ylabel('$x_2$')
ax.set_zlabel('$x_3$')
plt.show()
Programmons la méthode Monte Carlo :
def phi(x, c):
return (abs(np.prod(x)) < c)*1
def monte_carlo_estimate(n, d, c):
"""
n : taille de l'échantillon
d : dimension
c : le paramètre pour la région A
"""
rn = np.random.multivariate_normal(np.zeros(d), np.eye(d), n)
I_hat = 0
for i in range(n):
I_hat = I_hat + phi(rn[i,:], c)
I_hat = I_hat / n
return I_hat
def monte_carlo_estimate_list(n, d, c):
"""
n : taille de l'échantillon
d : dimension
c : le paramètre pour la région A
"""
rn = np.random.multivariate_normal(np.zeros(d), np.eye(d), n)
rv = [phi(rn[i,:], 0.9) for i in range(n)]
I_hat = np.cumsum(rv)
I_hat = I_hat/np.arange(1,len(rv)+1)
return I_hat
monte_carlo_estimate(n = 100, d = 3, c = 1)
monte_carlo_estimate(n = 100, d = 3, c = 0.5)
monte_carlo_estimate(n = 100, d = 3, c = 0.0001)
# lorsque c est petit, l'événement devient rare et l'estimation donne presque tout le temps 0
monte_carlo_estimate(n = 10000, d = 3, c = 0.0001)
# il est possible d'avoir une meilleure estimation en augmentant la taille de l'échantillon
N = 200
x = np.arange(1, N+1)
y = []
for n in x:
y.append( monte_carlo_estimate(n, d = 3, c = 0.5) )
plt.scatter(x,y)
plt.xlabel("n (taille de l'échantillon)")
plt.ylabel("$\hat I_n$")
plt.show()
x = np.arange(1, N+1)
y = []
for n in x:
y.append( monte_carlo_estimate(n, d = 3, c = 0.0001) )
plt.scatter(x,y)
plt.xlabel("n (taille de l'échantillon)")
plt.ylabel("$\hat I_n$")
plt.show()
Nous nous situons dans une situation où l'événement est rare et la méthode monte carlo "classique" n'est pas stable pour $n$ petit.
Idée : Prendre la statistique $h(X) = \sum\limits^{d}_{k=1}X^k$ comme variable de contrôle, pour deux raisons :
$h(X)$ est corrélé à $\{X \in A\}$
il est facile de calculer l'espérance de $h(X)$ : $\mu = \mathbb{E}[h(X)] = d\cdot \mathbb{E}[X^1] = d\cdot \sqrt{\frac{2}{\pi}}$
Jettons un coup d'oeil sur la forme de la fonction $\varphi$ :
def h(x):
return(sum(abs(x)))
def control_variates_estimate(n, d, c):
"""
n : taille de l'échantillon
d : dimension
c le paramètre pour la région A
"""
rn = np.random.multivariate_normal(np.zeros(d), np.eye(d), n)
fX = [phi(rn[i,], c) for i in range(n)]
fXbar = np.mean(fX)
hX = [h(rn[i,]) for i in range(n)]
mean_hX = d*np.sqrt(2/np.pi)
var_hX = np.var(hX)
# "[0][1]" car la fonction "cov" de numpy donne une matrice et l'élement que nous chercons se trouve dans la coord [0][1]
cov_hX_fX = np.cov(fX, hX)[0][1]
beta = - cov_hX_fX / var_hX
Z = fXbar + beta * np.mean(hX - mean_hX)
return(Z)
n = 200
d = 3
c = 0.01
Imc = []
Icv = []
for i in range(n):
Imc.append( monte_carlo_estimate(n, d, c) )
Icv.append( control_variates_estimate(n, d, c) )
plt.boxplot((Imc, Icv), labels = ["Monte Carlo", "Control Variates"])
plt.show()
Idée : Nous pouvons générer
Programmons tout d'abord la méthode Monte Carlo qui va nous servir de baseline
# monte carlo standard
def phi(x):
# x : un vecteur de dimension 2
return np.log(1 + abs(x[0] - x[1]))
def monte_carlo_estimate(n, d = 2, c = 0.5):
"""
n : taille de l'échantillon
d : dimension (2 par default)
c le paramètre pour la région A (0.5 par default)
"""
rn = np.random.multivariate_normal(np.zeros(d), np.eye(d), n)
index = (abs(np.prod(rn, axis=1)) < c)
rn = rn[index, :]
m = len(rn)
I_hat = 0
if m>0:
for i in range(m):
I_hat = I_hat + phi(rn[i,])
I_hat = I_hat / m
return(I_hat, m/n)
monte_carlo_estimate(100, d=2, c=0.01)
Proposition de simulation :
Simuler les $X_i \sim \mathcal{N}_d(0, I_d)$
Calculer la moyenne $\mu$ et la matrice covariance empirique $\hat \Sigma$ des $X_i$ appartenant à $A$
Faire de l'IS avec la loi auxiliaire $\mathcal{N}_d(\hat\mu, \hat \Sigma)$
Calculer $\hat\mu$ et $\hat \Sigma$ est une mauvaise idée car ces termes sont encore des estimateurs de MC standard sur des événements rares lorsque $c$ devient petit ...
Code :
from scipy import stats
def IS_gauss_estimate(n, d = 2, c = 0.5):
# etape 1
X = stats.multivariate_normal.rvs(np.zeros(d), np.eye(d), n)
index = (abs(np.prod(X, axis=1)) < c)
X = X[index, :]
m = len(X)
# etape 2
mu = np.mean(X, axis=0)
sig = np.cov(X.T)
# etape 3
Y = stats.multivariate_normal.rvs(mu, sig, n)
g = stats.multivariate_normal.pdf(Y, mu, sig)
f = stats.multivariate_normal.pdf(Y, np.zeros(d), np.eye(d))
index = (abs(np.prod(Y, axis=1)) > c)
f[index] = 0
m = sum(index)
W = f / g
phi_g = np.array(list(map(phi, Y)))
I_tilde = sum(phi_g * W) / sum(W)
return(I_tilde, m/n)
print(monte_carlo_estimate(10000, d=2, c=0.01))
print(IS_gauss_estimate(10000, d=2, c=0.01))
On change la loi auxiliaire / de proposition $g$ par la t-distribution qui est plus épaisse aux bords qu'une gaussienne
def IS_t_estimate(n, d=2, c=0.1, df=5):
rnd = np.random.RandomState()
# On genere une t de taille (n, d)
Xg = rnd.standard_t(df, size=(n, d))
# On calcule la densite de la loi t non normalisee
g = (1 + (Xg ** 2).sum(axis=1) / df) ** (- (2 + df) / 2)
# et celle de f (gaussienne tronquee)
f = np.exp(- 0.5 * (Xg ** 2).sum(axis=1))
index = (abs(np.prod(Xg, axis=1)) > c)
f[index] = 0.
weights = f / g
weights = np.nan_to_num(weights)
m = sum(index)
phi_g = np.array(list(map(phi, Xg)))
phibar = (phi_g * weights).sum() / weights.sum()
effective_n = weights.sum() ** 2 / (weights ** 2).sum()
return np.nan_to_num(phibar), m/n
n, d, c = 10000, 2, 0.001
print(monte_carlo_estimate(n, d, c))
print(IS_gauss_estimate(n, d, c))
print(IS_t_estimate(n, d, c))
n, d, c = 10000, 2, 0.0001
print(monte_carlo_estimate(n, d, c))
print(IS_gauss_estimate(n, d, c))
print(IS_t_estimate(n, d, c))
n, d, c = 1000, 2, 0.01
N = 50
data = np.zeros((N, 3))
data[:, 0] = np.array([monte_carlo_estimate(n, d, c)[0] for _ in range(N)])
data[:, 1] = np.array([IS_gauss_estimate(n, d, c)[0] for _ in range(N)])
data[:, 2] = np.array([IS_t_estimate(n, d, c)[0] for _ in range(N)])
plt.figure()
plt.boxplot(data, labels=["MC classique", "IS gaussian", "IS t"])
plt.show()
code : SOON ....